home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / a_utils / expanded.lha / test_suite / thex.C < prev   
C/C++ Source or Header  |  1992-03-19  |  16KB  |  515 lines
  1. //
  2. // Linear-Affine-Projective Geometry Package
  3. //
  4. // thex.C
  5. //
  6. // $Header$
  7. //
  8. // William J.R. Longabaugh 
  9. // University of Washington
  10. //
  11. // Thesis examples for the linear-affine-projective geometry
  12. // package described in William J.R. Longabaugh, "An Expanded
  13. // System for Coordinate-Free Geometric Programming", Master's 
  14. // thesis, University of Washington, 1992.
  15. //
  16. // Copyright (c) 1992, William J.R. Longabaugh
  17. //   Copying, use and development for non-commercial purposes permitted.
  18. //   All rights for commercial use reserved.
  19. //   This software is unsupported and without warranty; it is
  20. //   provided "as is".
  21. //
  22. // ***********************************************************************
  23.  
  24. #include "Lap.h"
  25. #include <math.h>
  26.  
  27. extern void thex1(void);
  28. extern void thex2(void);
  29. extern void thex3(void);
  30. extern void thex4(void);
  31. extern void thex5(void);
  32. extern void thex6(void);
  33. extern void thex7(void);
  34.  
  35. // ***********************************************************************
  36.  
  37. int main(void)
  38. {
  39.   thex1();
  40.   thex2();
  41.   thex3();
  42.   thex4();
  43.   thex5();
  44.   thex6();
  45.   thex7();
  46.   return (0);
  47. }
  48.  
  49. // ***********************************************************************
  50.  
  51. void thex1(void)
  52. {
  53. // Create the range and domain space for the map.  Note that
  54. // since this is a second-degree surface, the domain space
  55. // of the MAM is a cartesian product space with two components:
  56.  
  57.   ASpace param("Parameter space", 2, TRUE);
  58.   Simplex smp1 = param.StdSimplex();
  59.   ASpace domsp("Cart. prod. domain space",
  60.                SpaceList(param, 2), FALSE);
  61.  
  62.   ASpace rng("Range space", 3, TRUE);
  63.   Frame fr2 = rng.StdBasis();
  64.  
  65. // Construct the symmetry set specification.  Since this is
  66. // a fully symmetric map, the list has one entry (2):
  67.  
  68.   IntList symspec(1);
  69.   symspec[0] = 2;
  70.  
  71. // Construct the control net for the surface.  Note the ordering
  72. // of the points in the net.  List concatenation is used to
  73. // build a large net easily:
  74.  
  75.   APoint aa(fr2, ScalarList(0.0, 0.0, 0.0, 1.0));
  76.   APoint ab(fr2, ScalarList(1.0, 0.0, 0.7, 1.0));
  77.   APoint bb(fr2, ScalarList(2.0, 0.0, 0.0, 1.0));
  78.   APoint ac(fr2, ScalarList(0.5, 1.0, 0.7, 1.0));
  79.   APoint bc(fr2, ScalarList(1.5, 1.0, 0.7, 1.0));
  80.   APoint cc(fr2, ScalarList(1.0, 2.0, 0.0, 1.0));
  81.   
  82.   GeObList net = GeObList(aa, ab, bb) + GeObList(ac, bc, cc);
  83.  
  84. // Build the multi-affine map for the surface:
  85.  
  86.   MAM surfmap = MAM(domsp, symspec, BasisList(smp1), rng, net);
  87.  
  88. // Build a point in the parameter space for evaluating the map:
  89.  
  90.   APoint pt1(smp1, ScalarList(0.4, 0.3, 0.3));
  91.  
  92. // Evaluate the map to obtain a point on the surface.  Note
  93. // that the argument to the map is a tuple with the same
  94. // point repeated, so the result is on the surface:
  95.  
  96.   APoint argpt1(domsp, pt1, pt1);
  97.   APoint surfpt = surfmap(argpt1);
  98.  
  99. // Partially evaluate the map once on the point pt1 in the
  100. // parameter space.  The result is an affine map:
  101.  
  102.   GeObList partial(2);
  103.   partial[0] = pt1;
  104.   AffineMap tanmap = surfmap(param, partial);
  105.  
  106. // Mapping vectors in the parameter space through this map will
  107. // produce vectors tangent to the surface at the point surfpt.
  108. // Use vectors from the standard frame:
  109.  
  110.   AVector surftan1 = tanmap(param.StdBasis()[0]);
  111.   AVector surftan2 = tanmap(param.StdBasis()[1]);
  112.  
  113. // Take the cross product of the resulting vectors to obtain a
  114. // normal vector.  Normalize the length and cast it into a dual
  115. // vector.  Note that the argument to sqrt() must be cast to
  116. // a scalar manually; if the sqrt() function were not being used,
  117. // no casting would be needed:
  118.  
  119.   MLM rngcross = rng.GetSpace(TANGENT).CrossProduct();
  120.   MLM rnginner = rng.GetSpace(TANGENT).InnerProduct();
  121.   AVector nvec = rngcross(GeObList(surftan1, surftan2));
  122.   nvec = nvec / sqrt(rnginner(GeObList(nvec, nvec)).ToScalar());
  123.   Vector surfnorm = nvec.Dual();
  124.   cout << nvec;
  125.   cout << surfnorm;
  126.  
  127. // Use this to calculate the cosine of the angle of a unit
  128. // vector from the normal:
  129.  
  130.   AVector testvec(fr2, ScalarList(1.0, 0.0, 0.0, 0.0));
  131.   Scalar angle = surfnorm.Apply(testvec);
  132.  
  133. // Extract the coordinates of the surface point with respect to
  134. // the standard basis:
  135.  
  136.   ScalarList surfptCoords = fr2(surfpt);
  137.  
  138.   cout << "Cosine of angle: " << angle << "\n";
  139.   cout << surfptCoords;
  140. }
  141.  
  142. // ***********************************************************************
  143.  
  144. void thex2(void)
  145. {
  146. // Create the parameter spaces, range space, and domain space
  147. // for the map:
  148.  
  149.   ASpace paramX("Parameter space X", 1, TRUE);
  150.   Simplex smp1 = paramX.StdSimplex();
  151.  
  152.   ASpace paramY("Parameter space Y", 1, TRUE);
  153.   Simplex smp2 = paramY.StdSimplex();
  154.  
  155.   ASpace rng("Range space", 3, TRUE);
  156.   Frame fr1 = rng.StdBasis();
  157.  
  158. // The domain space of the MAM is a cartesian product space with
  159. // four components, two in each symmetry set:
  160.  
  161.   ASpace domsp("Cart. prod. domain space",
  162.                SpaceList(paramX, 2) + SpaceList(paramY, 2),
  163.                FALSE);
  164.  
  165. // Construct the symmetry set specification.  There are two
  166. // sets, each with two spaces (2,2):
  167.  
  168.   IntList symspec(2);
  169.   symspec[0] = 2;
  170.   symspec[1] = 2;
  171.  
  172. // Construct the control net for the surface.  Note the ordering
  173. // of the points in the net.  Use list concatenation to build
  174. // a large list easily:
  175.  
  176.   APoint aa_cc(fr1, ScalarList(0.0, 0.0, 0.0, 1.0));
  177.   APoint ab_cc(fr1, ScalarList(1.0, 0.0, 0.5, 1.0));
  178.   APoint bb_cc(fr1, ScalarList(2.0, 0.0, 0.0, 1.0));
  179.   APoint aa_cd(fr1, ScalarList(0.0, 1.0, 0.5, 1.0));
  180.   APoint ab_cd(fr1, ScalarList(1.0, 1.0, 1.0, 1.0));
  181.   APoint bb_cd(fr1, ScalarList(2.0, 1.0, 0.5, 1.0));
  182.   APoint aa_dd(fr1, ScalarList(0.0, 2.0, 0.0, 1.0));
  183.   APoint ab_dd(fr1, ScalarList(1.0, 2.0, 0.5, 1.0));
  184.   APoint bb_dd(fr1, ScalarList(2.0, 2.0, 0.0, 1.0));
  185.   
  186.   GeObList net = GeObList(aa_cc, aa_cd, aa_dd) +
  187.                  GeObList(ab_cc, ab_cd, ab_dd) +
  188.                  GeObList(bb_cc, bb_cd, bb_dd);
  189.  
  190. // Build the multi-affine map for the surface:
  191.  
  192.   MAM surfmap(domsp, symspec, BasisList(smp1, smp2), rng, net);
  193.  
  194. // Evaluate the map to obtain a point on the surface.
  195. // First build a point in each parameter space, then 
  196. // build a point tuple:
  197.  
  198.   APoint x(smp1, ScalarList(0.55, 0.45));
  199.   APoint y(smp2, ScalarList(0.3, 0.7));
  200.   APoint argpt1(domsp, GeObList(x, x) + GeObList(y, y));
  201.  
  202.   APoint surfpt = surfmap(argpt1);
  203.  
  204. // Extract the coordinates of the surface point with respect to
  205. // the standard basis:
  206.  
  207.   ScalarList surfptCoords = fr1(surfpt);
  208.  
  209.   cout << surfptCoords;
  210. }
  211.  
  212. // ***********************************************************************
  213.  
  214. void thex3(void)
  215. {
  216. // Represent a rational quadratic Bezier surface.  First create
  217. // the range and domain space for the map.  Both are
  218. // affine spaces:
  219.  
  220.   ASpace param("Parameter space", 2, TRUE);
  221.   Simplex smp1 = param.StdSimplex();
  222.   ASpace domsp("Cart. prod. domain space",
  223.                SpaceList(param, 2), FALSE);
  224.  
  225.   ASpace rng("Range space", 3, TRUE);
  226.   Frame fr2 = rng.StdBasis();
  227.  
  228. // Construct the symmetry set specification.  Since this is a 
  229. // fully symmetric map, the list has one entry (2):
  230.  
  231.   IntList symspec(1);
  232.   symspec[0] = 2;
  233.  
  234. // Construct the control net for the surface.  Since this is a
  235. // rational map, each control point has an associated
  236. // scalar weight, so build the corresponding list of scalars:
  237.  
  238.   APoint aa(fr2, ScalarList(0.0, 0.0, 0.0, 1.0));
  239.   APoint ab(fr2, ScalarList(1.0, 0.0, 0.7, 1.0));
  240.   APoint bb(fr2, ScalarList(2.0, 0.0, 0.0, 1.0));
  241.   APoint ac(fr2, ScalarList(0.5, 1.0, 0.7, 1.0));
  242.   APoint bc(fr2, ScalarList(1.5, 1.0, 0.7, 1.0));
  243.   APoint cc(fr2, ScalarList(1.0, 2.0, 0.0, 1.0));
  244.   
  245.   GeObList net = GeObList(aa, ab, bb) + GeObList(ac, bc, cc);
  246.   ScalarList weights = ScalarList(2.4, 1.3, -3.5) +
  247.                        ScalarList(5.1, 0.7, 4.6);
  248.  
  249. // Build the multi-affine map used to represent the surface.
  250. // The range of this map will be the linearization space of 
  251. // the affine space containing the surface.  The control net
  252. // of affine points is converted into a control net of vectors
  253. // in the linearization space: 
  254.  
  255.   VSpace tempRange = rng.GetSpace(LINEARIZATION);
  256.   GeObList vecnet(6);
  257.  
  258.   for (int i = 0; i < 6; i++) {
  259.     vecnet[i] = weights[i] * net[i];
  260.   }
  261.  
  262.   MAM surfmap(domsp, symspec,
  263.               BasisList(smp1), tempRange, vecnet);
  264.  
  265. // Build a point in the parameter space for evaluating the map:
  266.  
  267.   APoint pt1(smp1, ScalarList(0.2, 0.3, 0.5));
  268.  
  269. // Evaluate the map to obtain a point on the surface.  The 
  270. // MAM returns a vector in the linearization space.  This
  271. // is converted into a projective point in the projective
  272. // completion space (unless it is the zero vector).  The
  273. // projective point is then converted to an affine point
  274. // in the original range space unless it is a point at
  275. // infinity:
  276.  
  277.   APoint argpt1(domsp, pt1, pt1);
  278.   APoint surfpt = surfmap(argpt1).MapTo(VEC_EC)
  279.                                  .MapTo(PROJ_POINT)
  280.                                  .MapTo(AFF_POINT);
  281.  
  282. // Extract the coordinates of the surface point with respect to
  283. // the standard basis:
  284.  
  285.   ScalarList surfptCoords = fr2(surfpt);
  286.  
  287.   cout << surfptCoords;
  288. }
  289.  
  290. // ***********************************************************************
  291.  
  292. void thex4(void)
  293. {
  294. // Build the affine space and get a simplex and frame:
  295.  
  296.   ASpace workspace("Working space", 3, TRUE);
  297.   Simplex smp1 = workspace.StdSimplex();
  298.   Frame fra1 = workspace.StdBasis();
  299.  
  300. // Define the face with three points, and three vectors
  301. // representing normals at each vertex.  Convert each normal
  302. // vector into its corresponding dual vector (linear
  303. // functional):
  304.  
  305.   APoint pt1(fra1, ScalarList(7.2, 1.4, 5.5, 1.0));
  306.   APoint pt2(fra1, ScalarList(3.2, 5.4, 9.8, 1.0));
  307.   APoint pt3(fra1, ScalarList(1.2, 4.1, 6.7, 1.0));
  308.   AVector nv1(fra1, ScalarList(2.5, 3.6, 1.8, 0.0));
  309.   AVector nv2(fra1, ScalarList(6.6, 3.4, 1.6, 0.0));
  310.   AVector nv3(fra1, ScalarList(7.2, 5.1, 3.9, 0.0));
  311.   Vector norm1 = nv1.Dual();
  312.   Vector norm2 = nv2.Dual();
  313.   Vector norm3 = nv3.Dual();
  314.  
  315. // We wish to calculate the linear interpolation of the normal
  316. // for any given point on the face.  This is an affine map from
  317. // 2D affine subset defined by the face to the affine subset
  318. // defined by the three vectors in the dual space.  Define the
  319. // domain and range subsets.  Remember that these last the life
  320. // of the program (their storage space is not released):
  321.  
  322.   ASubSet face("face", workspace, GeObList(pt1, pt2, pt3));
  323.   ASubSet norms("norms", workspace.GetSpace(TANG_DUAL),
  324.                 GeObList(norm1, norm2, norm3));
  325.  
  326. // Define the map that does interpolation of the normals:
  327.  
  328.   AffineMap interp(face, GeObList(pt1, pt2, pt3),
  329.                    norms, GeObList(norm1, norm2, norm3));
  330.  
  331. // Create a point on the face and calculate the corresponding
  332. // normal:
  333.  
  334.   APoint testpt = .6 * pt1 + .2 * pt2 + .2 * pt3;
  335.   Vector intnorm = interp(testpt);
  336.  
  337.   cout << intnorm;
  338.  
  339. // Note that the interpolation map is invertible.  Given a
  340. // normal from the affine subset, get the corresponding
  341. // point on the face:
  342.  
  343.   AffineMap invint = interp.Inv();
  344.   Vector testnorm = .5 * norm1 + .1 * norm2 + .4 * norm3;
  345.   APoint intpt = invint(testnorm);
  346.  
  347.   cout << intpt;
  348. }
  349.  
  350. // ***********************************************************************
  351.  
  352. void thex5(void)
  353. {
  354. // Create an affine space.  The domain of the determinant is the
  355. // tangent space, while the range is the real numbers:
  356.  
  357.   ASpace workspace("Working space", 3, TRUE);
  358.   Simplex smp1 = workspace.StdSimplex();
  359.   Frame fra1 = workspace.StdBasis();
  360.   VSpace wsTangent = workspace.GetSpace(TANGENT);
  361.   VBasis bas1 = wsTangent.StdBasis();
  362.  
  363. // Construct the symmetry set specification.  Determinants in
  364. // 3D spaces take 3 arguments, and the negative sign indicates
  365. // the map is antisymmetric:
  366.  
  367.   IntList symspec(1);
  368.   symspec[0] = -3;
  369.  
  370. // Define the MLM.  The ``control net'' for the map is a single
  371. // vector: the scalar 1.0 in the reals.
  372.  
  373.   VSpace domsp("Cart. prod. domain space",
  374.                SpaceList(wsTangent, 3), FALSE);
  375.   MLM det(domsp, symspec, BasisList(bas1),
  376.           Reals, GeObList(Vector(1.0))); 
  377.  
  378. // Build a volume out of vectors:
  379.  
  380.   GeObList volume(3);
  381.  
  382.   APoint pt1(fra1, ScalarList(7.2, 1.4, 5.5, 1.0));
  383.   APoint pt2(fra1, ScalarList(3.2, 5.4, 9.8, 1.0));
  384.   APoint pt3(fra1, ScalarList(1.2, 4.1, 6.7, 1.0));
  385.   APoint pt4(fra1, ScalarList(3.4, 9.5, 7.1, 1.0));
  386.  
  387.   volume[0] = pt2 - pt1;
  388.   volume[1] = pt3 - pt1;
  389.   volume[2] = pt4 - pt1;
  390.  
  391. // Find the signed volume represented by the vectors.  The result
  392. // of the map is a vector in the special space "Reals", which
  393. // is manually cast to a scalar:
  394.  
  395.   Scalar res = det(volume).ToScalar();
  396.  
  397.   cout << "Determinant: " << res << "\n";
  398. }
  399.  
  400. // ***********************************************************************
  401.  
  402. void thex6(void)
  403. {
  404. // Create a 3D affine space:
  405.  
  406.   ASpace workspace("Working space", 3, TRUE);
  407.   Simplex smp1 = workspace.StdSimplex();
  408.   Frame fra1 = workspace.StdBasis();
  409.  
  410. // Get a normal and a point for the hyperplane:
  411.  
  412.   AVector nv(fra1, ScalarList(1.0, 1.0, 1.0, 0.0));
  413.   Vector norm = nv.Dual();
  414.   APoint pt(fra1, ScalarList(0.5, 0.5, 0.5, 1.0));
  415.  
  416. // Define the map by describing the action of the hyperplane
  417. // on the standard simplex:
  418.  
  419.   GeObList res(workspace.Dim() + 1);
  420.   for (int i = 0; i <= workspace.Dim(); i++) {
  421.     res[i] = norm.Apply(smp1[i] - pt);
  422.   }
  423.  
  424.   AffineMap hyp(smp1, Reals.FullSet(), res);
  425.  
  426. // Given a line segment that intersects the hyperplane, compute
  427. // the point of intersection:
  428.  
  429.   APoint endpt1(fra1, ScalarList(-2.3, -7.8, -11.2, 1.0));
  430.   APoint endpt2(fra1, ScalarList(4.6, 5.8, 6.7, 1.0));
  431.   Scalar d1 = fabs(hyp(endpt1).ToScalar());
  432.   Scalar d2 = fabs(hyp(endpt2).ToScalar());
  433.  
  434.   APoint intsct = (d1 / (d1 + d2)) * endpt2 +
  435.                   (d2 / (d1 + d2)) * endpt1;
  436.  
  437. // This point is in the hyperplane:
  438.  
  439.   MLM cross = workspace.GetSpace(TANGENT).CrossProduct();
  440.   Scalar outp =
  441.        hyp(pt + cross(GeObList(nv, fra1[0]))).ToScalar();
  442.  
  443.   cout << "Output: " << outp << "\n";
  444.  
  445.   cout << intsct;
  446.   outp = hyp(intsct).ToScalar();
  447.   cout << "Output: " << outp << "\n";
  448. }
  449.  
  450. // ***********************************************************************
  451.  
  452. void thex7(void)
  453. {
  454. // The map carries projective points from a 2-dimensional space
  455. // to a 1-dimensional subset of another 2-dimensional space:
  456.  
  457.   PSpace domsp("Domain space", 2);
  458.   HFrame fra1 = domsp.StdBasis();
  459.   PSpace rngsp("Range space", 2);
  460.   HFrame fra2 = rngsp.StdBasis();
  461.  
  462. // Construct the domain subset, which can be viewed as modeling
  463. // a pencil of lines through the point p1.  Points p1, p4, and
  464. // p2 define a plane in the embedding space.  The point p1 is
  465. // then removed from this plane, creating the subset:
  466.  
  467.   PPoint p1(fra1, ScalarList(6.7, 2.3, 7.5));
  468.   PPoint p2(fra1, ScalarList(3.6, -1.4, 2.2));
  469.   PPoint p3(fra1, ScalarList(1.1, 0.2, 1.0));
  470.   PPoint p4(fra1, ScalarList(5.6, 3.4, -9.1));
  471.   PSubSet domsub("Domain subset", domsp,
  472.                  GeObList(p1), GeObList(p4, p2));
  473.  
  474. // Create the range subset, which consists of all the points on
  475. // the line pr1pr2:
  476.  
  477.   PPoint pr1(fra2, ScalarList(10.7, 3.4, 1.1));
  478.   PPoint pr2(fra2, ScalarList(5.6, -1.4, 12.2));
  479.   PPoint pr3(fra2, ScalarList(16.3, 2.0, 13.3));
  480.   PSubSet rngsub("Range subset", rngsp, GeObList(pr1, pr2));
  481.  
  482. // Create the non-invertible map:
  483.  
  484.   ProjectiveMap pm1(domsub, GeObList(p4, p2, p3),
  485.                     rngsub, GeObList(pr1, pr2, pr3));
  486.  
  487. // Send a point through the projective map.  The algebraic
  488. // operations on projective points are performed in the
  489. // context of the standard affine subset defined on the
  490. // projective space:
  491.  
  492.   PPoint ppo = pm1((p1 + p2) / 2.0);
  493.  
  494. // Obtain the homogeneous coordinates of the resulting point:
  495.  
  496.   ScalarList ppoCoords = fra2(ppo);
  497.   cout << ppoCoords;
  498.  
  499.   ppo = pm1(p4);
  500.   ppoCoords = fra2(ppo);
  501.   cout << ppoCoords;
  502.  
  503.   ppo = pm1(p2);
  504.   ppoCoords = fra2(ppo);
  505.   cout << ppoCoords;
  506.  
  507.   ppo = pm1(p3);
  508.   ppoCoords = fra2(ppo);
  509.   cout << ppoCoords;
  510.  
  511.   ppo = pm1((p1 + p4) / 2.0);
  512.   ppoCoords = fra2(ppo);
  513.   cout << ppoCoords;
  514. }
  515.